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On the basis of a gaussian quasi-chemical model of hydration, a model of non van der Waals 
character, we explore the role of attractive methane-water interactions in the hydration of methane 
and in the potential of mean force between two methane molecules in water. We find that the hy- 
dration of methane is dominated by packing and a mean-field energetic contribution. Contributions 
beyond the mean-field term are unimportant in the hydration phenomena for a hydrophobic solute 
such as methane. Attractive solute-water interactions make a net repulsive contribution to these 
pair potentials of mean force. With no conditioning, the observed distributions of binding energies 
are super-gaussian and can be effectively modeled by a Gumbel (extreme value) distribution. This 
further supports the view that the characteristic form of the unconditioned distribution in the high-e 
tail is due to energetic interactions with a small number of molecules. Generalized extreme value 
distributions also effectively model the results with minimal conditioning, but in those cases the 
distributions are sufficiently narrow that details of their shape aren't significant. 



I. INTRODUCTION 

It is generally accepted that hydrophobic effects are 
an essential contribution to the stability of the func- 
tional structures of soluble proteins!^ After Frank and 
Evangel and then Kauzmannpl hydrophobic effects have 
been identified by distinctive unfavorable entropies of hy- 
dration and large hydration heat capacities. That has 
made development of a fully defensible molecular theory 
of hydrophobic effects an outstanding challenge for sta- 
tistical thermodynamics of solutions. 

Over the past decade, however, such a defensible 
molecular t heory has begun to take shapeP^ Those 
development j 6 * 7 * 8 -^ started principally from the scaled- 
particle models but also de fined a connection to the 
Pratt-Chandler theory. El 12 l 13 l In view of those develop- 
ments, consolidation of previous work is to be expectedP 
This paper addresses the contributions of attractive in- 
teractions, e.g., dispersion interactions, on the potentials 
of the average forces between small hydrophobic solutes 
in water. That is a topic of specific interest that has been 
left unresolved. 

Much of that specific interest is associated with the 
Pratt-Chandler theory which was the first serious the- 
oretical prediction of the potential of mean force, p mf. 
between simple hydrophobic solutes in waterPA paper^l 
which followed the initial Pratt-Chandler predictions 
treated the role of attractive solute-solvent interactions, 
and addressed an apparent disagreement with experi- 
ment that had surfaced.^ By now, a significant body 
of simulation data has accumulated on those pmfs. In 
the following, we briefly note those results salient to our 
goals. But even more, a recent theoretical insightPSI seems 
to provide just the analysis tool needed to organize and 
clarify the accumulated simulation data, and to clarify 
the theoretical understanding of these pmfs. 



The apparent disagreement involved the temperature 
dependence of the osmotic second virial coefficient 

B2 = -\J {e- w ^/ kT - l) 47rr 2 dr , (1) 

where waa (t) is the pmf between two hydrophobic so- 
lutes A, assumed to be spherical in the original Pratt- 
Chandler theory which treated hard-sphere solutes. The 
sign of B 2 predicted for hard-sphere solutes by the Pratt- 
Chandler theory depended on the size of the solute, the 
predictions for larger solutes was attractive B 2 <0. But 
the magnitude of B 2 decreased with increasing temper- 
ature, so the effective attractive interactions were pre- 
dicted to become less attractive as the temperature was 
raised. The general expectation for hydrophobic inter- 
actions always had been that B 2 should be negative, 
indicating a net attraction, and increasing in strength 
with increasing temperature. Subseque nt ex periments 
for the non-spherical A = GqHq (benzene) ! 17 * 18 ! and for A 
= CgF 6 (perfluoro-benzene^PIl were consistent with the 
general expectation for hydrophobic interactions: this is 
a favorable and endothermic association, i.e., this is an 
entropy driven association. Though later results have 
revised that endothcrmicity for benzene association to 
smaller valuesj22l simulation results on model spherical 
solutes in water agree with the prior expectation for hy- 
drophobic interactions ! 21 * 22 ^ 

We note also that study of the benzene-benzene pmf 
by direct numerical simulation has pro duced res ults in 
rough agreement with the experiments! 23 * 24 * 25 * 26 ! These 
are sufficiently challenging calculations, however, that 
those simulation results may not be definitive for B 2 . 
A significant result of those calculations is that the 
T-shaped benzene-benzene c ontacts are more favorable 
than stacked configurations.^ 25 * 26 * This seems to be 
qualitatively different from the case studied with purely 
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repulsive model plates stacked.^ This detail suggests a 
general conclusion that attractive solute-solvent interac- 
tions can play a qualitative role in establishing the most 
probable contacts. That realistically modelled benzene 
dimers and toluene dimers may differ significantly in this 
regard^ indicates yet again that realistic details of the 
cases considered may decisively affect these contacts . 

For spherical hydrophobic solutes, the case of Kr(aq) 
has received special attention. Whether B 2 > fol- 
lows from the experimental solubilities of Kr(aq) is 
undecidedp^l but relying on simulation results again, we 
do know that it is a realistic possibility.^ Those simula- 
tion results, furthermore, do suggest that the the sign of 
B 2 can change with the strength of solute-solvent attrac- 
tive interactions. With the potential of mean force for a 
C60 pairj^fll in a different size regime, attractive interac- 
tions make a dominating contribution that is repulsive in 
character. 

The work of Pangali et aZ.pH indicated that the original 
Pratt-Chandler theory was valid for waa ( r ) in the case 
where the predicted contribution of attractive interac- 
tions was small J 1 ^ Furthermore, the theory of Ref. 14 was 
observed by Smith & Haymell 3 - 2 -' on the basis of molecular 
simulation to be strikingly accurate for waa ( r ) point- 
wise in the Lennard- Jones model case there considered 
for which the predicted contribution of attractive inter- 
actions was large. For that case B 2 > Op^l Nevertheless, 
the underlying Pratt-Chandler theory was the source of 
an inaccurate association of that Lennard- Jones model 
with CH4(aq). As an example of the type of inaccuracy 
that should be expected, the original Pratt-Chandler the- 
ory is approximate even at small solute sizes for the case 
of a hard-sphere sol ventPS Thus the accuracy of the pmf 
result for that case 32 was not developed further. 

The observations above suggest the following general 
hypotheses: (a) the balance of attractive interactions can 
play a significant role in hydrophobic interactions; (b) 
with attractive interactions included the original Pratt- 
Chandler theory can be accurate for the hydrophobic 
interactions expressed by waa(t). The original Pratt- 
Chandler theory for hydrophobic hydration problems was 
less accurate. But the modern revisions of that theory ad- 
dress those hydration problems first and largely on the 
basis of scaled -par ticle approaches^ which are by now 
well-developed. 1 ™ ( c ) B 2 > is a realizable possibility, 
but probably not for the cases of general experimental 
interest. 

In response to the apparent experimental disagreement 
noted above, a variety of complications were considered: 
(a) context hydrophobicity which, for example, might be 
expressed in experiments on solutes such as methanol in 
aqueous solution.^ (b) non-spherical shapes of molecules 
of experimental interest,^ and (c) the influence of attrac- 
tive solute-solvent attractive interactions.^ These are all 
valid general concerns. But the post hoc conclusion has 
been that the underlying Pratt-Chandler theory is suf- 
ficiently unconvincing that extensions and elaborations 
have not been compelling. 



Though the intuitive proposal of Ref. 14 for inclu- 
sion of attractive intera ctions worked-out accurately in 
some cases testedj^^l there is a logical complication. 
It is common to assume that attractive interactions 
may be included perturbatively after packing problems 
dominated by repulsive solute-solvent interactions are 
solved. Then it is expected that attractive interac- 
tions may change thermodynamic properties significantly 
while making small alterations to structural characteris- 
tics. Being an integrated quantity, B 2 is a thermody- 
namic property and should be expected to be sensitive 
to inclusion of attractive interactions, from this perspec- 
tive. But the assumed insensitivity of structural charac- 
teristics has not been so carefully examined; a physical 
comment about the complicated role of hard-core models 
of hydrophobic solutes was attempted in Ref. 5, and fur- 
ther observations on the role of attractive interactions in 
hydrophobic solubility models were made by Paschek.^ 

A first perspective on the complications that can 
arise in starting from hard-core model solutes is that 
those complications are precursors of the drying that 
is expected on large length scales. ^ Dr ying has been 
discussed broadly in recent timespSMl Similarly, the 
startling success of the information the ory approach for 
hard-sphere solutes of intermediate siz d 8 * 40 ! derives from 
a balance of errors, approximate treatment of packing 
problems balanced by neglect of incipient drying. 4 An ad- 
ditional perspective is that the difficulties of the original 
Pratt-Chandler theory are associated with detailed prob- 
lems of Percus-Yevick-analogue theories as was noted 
aboveP^ Multiple complications of different types makes 
correcting a reference (hard-core) theory tricky. Never- 
theless, the status of the Pratt-Chandler theory for hard- 
sphere solutes in water has changed^, and the amended 
theory is now a compelling approximate theory — a 
quasi-chemical theory^ — relying o n and acknowledging 
copious empirical input pi ° l 41 l 42 l 43 l 

With this context, the work of Ref. 16 provides a key 
to theoretical assessment of the role of attractive interac- 
tions for these problems. That theory utilizes a result for 
the hard-core problem to address cases with attractive in- 
teractions, but the logical status of the hard-core result is 
not that of a reference system. There is no uncontrolled 
assumption of linear perturbative treatment of attractive 
interactions. Indeed, below we calculate non-linear con- 
tributions of two distinct types beyond non-hard-core in- 
teractions and an intuitively recognizable mean-field con- 
tribution, as is discussed below. Characterization and 
testing of this alternative theory for such a pmf is thus 
the target of this paper. 



II. THEORY 

We recapitulate the gaussian model of hydration and 
in the process also set the notation. Then we consider 
the extension of the theory to describe the potential of 
mean force between two methane molecules in water. 
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On the basis of the inverse form of the potential 
distribution theorem, the excess chemical potential of 
methane(aq), /i^J, is given by 



<e*> 



P M (e)e^ de 



(2) 



(3 = l/ksT, where T is the temperature. (. . .) specifies 
averaging on the basis of the probability density function 
-Pm(e) = — AUm)), which defines the distribution 
of binding energies of the solute with the solvent. For 
a particular configuration, the binding energy AUm — 
Un+i — Un — Ui is the difference in the potential energy of 
the whole system (f/jv+i) and the sum of the (decoupled) 
potential energies of the solution (Un) and the solute 
(Ui). Eq. ([2| is independent of the simulation ensemble 
in the macroscopic limit and an excellent approximation 
for system sizes considered in this worlP^. 

Notice that the exponential weighting exhibited by 
Eq. (12| emphasizes high binding-energy contributions to 
Pm(£)- These high energy contributions reflect collisions 
at short-range between the solute and the solvent parti- 
cles. These high-e (low-probability) features also render 
the direct application of Eq. Q problematic. 

Instead we regularize of the statistical problem as fol- 
lows: Consider a hard-core (HC) solute that excludes 
water oxygen atoms from a spherical region of radius A 
centered on the methane. A*hC' ^ ne excess chemical po- 
tential of this hard-core solute, is assumed to be known. 
On the basis of simulation data obtained from a calcula- 
tion on methane-water system with realistic interactions, 

Mhc is s iven b y 



-/3(£hc-s)\ 



(3) 



£hc is the binding energy of the hard-core model solute; 
this is zero when no water molecules overlap the hard- 
core model solute and is infinite otherwise, e is the bind- 
ing energy of M. The use of Eq. (|} requires eliminating 
configurations in which a solvent molecule penetrates the 
assigned hard-core exclusion volume in a simulation of M 
in water. The fraction of configurations thus selected is 
Pm(«a = 0). Then 

/?Mm = Mo + lnp M K = 0) 

+ ln /p M (£|nA = 0)e /3£ de . (4) 



Here Pm(£|"-a = 0) is the binding energy distribution con- 
ditional on the event that there are no overlaps with the 
defined hard-core model solute. The virtue of the reg- 
ularization is that the troublesome high-e features are 
accounted for by /i£f c and %>yi{ri\ — 0). The process is 
one of pushing the water molecules away from the solute. 
Thus by controlling the spatial distribution of molecules 
we temper the binding energies. The binding energies 
remaining after rcgularization are less energetic. They 
are composed of contributions from numerous sources, 
distant and thus weakly correlated. It is expected that 



Pm(e|^a = 0) will be well described by a normal dis- 
tribution. If Pyii^nx = 0) is modelled as a Gaussian 
distribution, then 

/?Mm « /Wc + lnp M («A = 0) 

+ P(e\n x = 0) + ( ^(Se 2 \n x = 0) . (5) 

This approach has been tested on the case of hydrophobic 
CF4(aq55_and on the case of liquid water in more than 
one way! 45 * 46 ] The balance between the packing and chem- 
ical contributions to Eq. ([5]), /i£f c and kBT\npM(n\ = 0), 
respectively, and particularly with changes in the volume 
of the defined inner-shell has been considered explicitly. 
A detailed point of note is that if the inner-shell is cho- 
sen small enough, then lnpyi(n\ = 0) s» 0, but the pack- 
ing contribution is non-zero. Then as the volume 
of the inner-shell is reduced, the packing contribution 
should be taken to be the largest value consistent with 
lnpM( n A = 0) ~ to sufficient accuracy. Thus, for exam- 
ple, with this prescription the gaussian estimate of the 
free energy of K + (aq) 47 can be slightly but distinctly im- 
proved by adding the appropriate packing contribution. 

With the free energy Eq. (J5j) in hand, the entropy can 
be accessed through the thermodynamic identity 



T 



PS 
dn M 



T,p,nw 



d(E) 
dn M 



T,p,Tlw 

d(V) 
dn M 



Mm 



(6) 



The ideal contributions to these thermodynamic quan- 
tities aren't specifically interesting, so we subtract the 
corresponding relation that obtains when the interaction 
potential energies vanish. Recalling that 



d(V) 
dn M 



- ideal 



P 



(7) 



yields 



T 



dS ex 
dn M 



d(U) 
dn M 

d(V) 
dn M 



Mm 



(8) 



Here (U) is the expected value of the potential energy 
of the system. The contribution of the partial molar en- 
ergy involves the mean binding energy — without the 
conditioning discussed above — but also the alteration 
of solvent-solvent interactions due to the presence of the 
solute. The significance of this term can be appreciated 
by considering the case of hard-sphere model methane. 
Then the theory Eq. (|5| without the outer-shell interac- 
tion contributions is transparently correct, and the par- 
tial molar energy contribution to Eq. (|8| involves only 
changes in solvent-solvent interactions. 
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Similarly, it is helpful to discuss the contribution to 
Eq. ^ associated with the partial molar volume. The 
magnitudes of solute partial molar volumes for hydropho- 
bic solutes are similar to, typically somewhat smaller 
than, molecular van der Waals volumes. For moderate 
pressures and in view of the the typical low solvent com- 
pressibility, the contribution from the actual partial mo- 
lar volume is negligible. The net contribution from the 
excess partial molar volume in Eq. ^ amounts invari- 
ably to —k^T, the specific value of the full partial molar 
volume in Eq. (|8| being numerically irrelevant. 



III. METHANE-METHANE POTENTIAL OF 
MEAN FORCE 



In this case, the second M molecule is physically present 
at the radius r in carrying-out the average and thus the 
subscript zero (0) is absent. The numerator in Eq. (12 1 
is the average for a solitary M molecule, or for the second 
M molecule infinitely removed from the first one located 
at the origin. The desired pmf is then 

- kTlnguu (r) = wmmW + AwmmW, ( 13 ) 

which defines AwmmM, an( l then further denotes 

AwMM(r) + ^ = kT In (e? AU ™\r) . (14) 

Following the steps that lead to Eq. ^ , here we condi- 
tion simulation data for the case of two actual M present 
and obtain 



Now let us consider how to use this statistical thermo- 
dynamic model to analyze the methane-methane pmf. 
Consider a first methane molecule centered at the origin. 
In fact, we will use a united-atom representation of the 
methane molecule. Then focus on the distribution of a 
second methane molecule in the field of the first. That 
distribution can be analyzed beginning from the princi- 
ple of constancy of the chemical potential in non-uniform 
system: 
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Pu(r) 



e -/3uMM(r) ^ e -/3Al/ M 



/ / -/3AU M \ r \ " 

-Pu MU (r) I \ r /a 



-0AU*, 



(9) 



Here mmmW is the assumed M-M pair interaction po- 
tential energy function, and Am is the thermal deBroglic 
wavelength for a methane molecule. The brackets (. . .) 
indicate a test particle average, i.e. averaging in the ab- 
sence of coupling between the medium and the second 
M molecule. AUm is the binding energy of the second 
M molecule with the rest of the solvent medium. The 
conditional test-particle average (. . . \r) is for the sec- 
ond M at position r. The unconditional average will be 
practically equal to the conditional average for positions 
r well separated from the first M molecule at the origin. 
The replacement 

for radial characteristics of the density profile provides 
the desired radial distribution function jmm (r) according 
to 



// e -/JAtf M | r V 

o MM (r) = e""° MuW - — 




(11) 



The inverse forrrPH corresponding to Eq. (Ill is 



5mm (r) = e-^ MM « 



J3AUi 



) 



/ e/ 3AC/ M | r ) 



(12) 



AwmmM + Mm ~ MHc( r ) + k BTlnp M (n x = 0\r) 

+ (e\r,n\ = 0) + (6e 2 \r,n x = 0) /2kT . (15) 

(The averaging is over a probability density distribution 
-Pm^I?") = (S(e — AUu)\r).) These formulae apply to the 
second M molecule physically present at radius r with 
the first M molecule located at the origin, and the indi- 
cated binding energies e are for the second M molecule. 
Nevertheless, 



^ x c (r) = -k B Tlnp(n x = 0\r) 



(16) 



interrogates spontaneous molecular scale cavity forma- 
tion in the neighborhood of a single M molecule, the first 
M molecule which is located at the origin. 

Eq. (15 1 is clearly approximate, but, as in Eq. |J the 
binding energy distribution of only one M molecule is 
regularized. Further it is a physically motivated theory 
and also systematically organized. Thus, it has the ad- 
vantages of physical theories, e.g. providing physical ex- 
planations of the observed behaviors. More specifically it 
is directly associated with the physical theory Eq. ^ , so 
this approach also offers an internal calibration, in con- 
trast to calculations of mean forces that determine the 
spatially varying pmfs but not a constant of integration. 



IV. METHODS 

We utilize computational results from simulations in- 
volving zero, one, and two methane molecules. A united 
atom representation of methane (Lennard- Jones e = 
0.294 kcal/mole and a — 3.73 A) was adopted. The 
SPC/E water model was used in all the simulations with 
the rigid H2O structure of each water molecule controlled 
by the SHAKE method. 

Classical molecular dynamics simulations were per- 
formed with the NAMD 2.6 code with constant NPT pro- 
cedures. The pressure was held constant at 1 bar using 
the Langevin piston method with the Langevin piston 
period set at 200 fs with a decay constant of 100 fs. The 
temperature of 298 K was controlled by Langevin dy- 
namics applied to the oxygen atoms. The temperature 
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of the Langevin piston was set equal to the temperature 
of the Langevin thermostat. The damping coefficient of 
the thermostat was set at 1 ps" 1 . 

A simulation of pure water was conducted for a system 
of 512 water molecules. The Lcnnard- Jones interaction 
was terminated at 10.5 A by smoothly switching to zero 
starting at 10 A. Electrostatic interactions were treated 
with the particle mesh Ewald method with a grid spac- 
ing of approximately 0.75 A. A well-equilibrated box of 
water molecule was used and a further equilibration for 
250 ps was conducted. The production run lasted 750 ps 
and configurations were saved every 0.2 ps for analysis. 
/3/iHc was evaluated by trial insertions into liquid wa- 
ter. Similarly we estimated the test particle distribution 
P^P(s) — (8{e — AU))o on the basis of about 2.7 million 
trial insertions of the methane into these water systems. 

For the methane in water simulation, a methane 
molecule was placed in the center of the water box. For 
convenience in analysis, the methane was held fixed. The 
system was extensively equilibrated for over 2 ns. The 
production run lasted 4 ns and configurations were saved 
every 0.5 ps. FIG. [T] shows the CO radial distribution 
function observed. 
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FIG. 1: The radial distribution function of water oxygen 
atoms around methane. 3.3 A, the smallest distance for which 
SCO (r) = 1, is chosen as the radius A of the inner shell. 

Pm(z) wa s assessed and its overlap with Pj^\e) was 
used to estimate according to the standard relation 
fcrjT 4n P M (e) / Pm\e) — — e + /i^. A conditioning ra- 
dius of 3.3 A was adopted in considering Pm(e|^a = 0). 
Thus only configurations in which no water molecules 
come closer than 3.3 A to the methane are retained for 
analysis. The same conditioning radius is used in all the 
gaussian models below. FIG. [2] shows the unconditional 
and conditional energy distributions. 

To compute n^ c (r) we adopted the following ap- 
proach. For each radial distance from the methane, 572 
points were placed to cover a sphere uniformly. Observa- 
tion spheres of radius 3.3 A were placed at each point and 
the water occupancy statistic compiled. Statistics were 
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FIG. 2: Distribution of binding energies of methane to water. 
The open circles are the unconditional distribution Pm(e). 
The filled circles are the conditional distribution Pm(e|^a = 0) 
obtained from the sample in which no water O-atoms occupy 
the volume defined by distances < 3.3 A from the methane 
center. The dashed curve is a gaussian model for the condi- 
tional distribution. The solid curves are generalized extreme 
valutP distributions fitted to the data. For the conditioned 
distribution this is the Frechet extreme value distribution 
lnP( £ ) = (n-l)ln(l-(^))-(l-(^)) n -ln6with 
n = 6.34, a = -3.76, and b = 0.37. The model plotted for 
the unconditioned distribution is the Gumbel extreme value 
distribution discussed in the text. 



accumulated from 8000 frames for about 4.6 million sam- 
ple points. p(0|r), the probability that the observation 
volume is empty immediately gives /iHc( r ) according to 
Eq. (pi 



The theory above was tested against numerically exact 
results obtained by an overlap method as follows. Firstly, 
we obtained the distribution Pm (e|r) of binding energies 
of one of two methane molecules that are present in a 
water system, the two methane molecules being placed 
along the box diagonal at various separations, r, ranging 
from 4 A to 8 A in Ar = 0.2 A intervals. The molecules 
were held fixed and the system equilibrated for over 1 ns. 
A production run 750 ps long was conducted with con- 
figurations saved every 0.2 ps. Secondly, we obtained the 
distribution pj^ (e\r) of a test molecule positioned at r 
in a corresponding simulation of water and one methane 
molecule. These distributions are related by 



-Pm (e|r) 
Wr) 



_ g-(e-Aw MM (r)-/i|J)//sT 



(17) 



Thus, evaluation of these distributions accurately enough 
to determine this ratio through an intermediate e regime 
determines Awmm(^) + Mm- 
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V. RESULTS AND DISCUSSION 

FIG. [T] shows the distribution of water oxygen atoms 
radially from the methane carbon. A = 3.3 A is the small- 
est distance for which gco (A) = 1. Below this radius, the 
influence of repulsive interaction dominates the interac- 
tions, and the conditional density of water becomes low. 
We define the radius of the hard sphere to be this value, 
and the probability pu{n\ = 0) follows from that. Here 
that probability is about PM( n \ = 0) ~ 0.6, so the con- 
ditioning eliminates only about half of the total number 
of configurations. The mean number of water O- atoms 
with r < A is also about 0.6. From this we can conclude 
that n\ = 1, and n\ = 2 are about equally likely with 
probability of 0.2. 

This minimal conditioning, nevertheless, has a dra- 
matic affect on the distribution of binding energies as 
FIG.EJshows. The extended tail seen there in the uncon- 
ditional distribution reflects the high-energy short-range 
collisions between the methane solute and the nearest 
few water molecules, as discussed above. The condi- 
tioning excludes those high-e interactions and the re- 
maining interactions conform to the anticipated gaussian 
distribution. For the conditional distribution (e\n\ — 
0) = —3.6 kcal/mole, whereas (e) = —3.2 kcal/mole. 
Based on the fraction of configurations that contribute 
to the conditional distribution, evacuation of the inner- 
shell requires the modest energy — fceTlnpM^A = 0) = 
0.29 kcal/mole. The latter quantity is the chemical con- 
tribution identified in earlier quasi-chemical generaliza- 
tions of the potential distribution theorem! 16 * 45 ! The fluc- 
tuation contribution (5e 2 \n\ = 0) /2kT = 0.13 kcal/mole. 
This is about 4% of the mean binding energy. Thus 
a mean-field approximation is adequate to describe the 
hydration of methane. By trial insertions, we estimate 
Mhc = 6.08 kcal/mole, in good agreement with the 
value of 6.2 kcal/mole predicted by revised scaled par- 
ticle theory^. Thus on the basis of Eq. (i^ = 
2.32 kcal/mole, in excellent agreement with the value 
of 2.38 kcal/mole obtained on the basis of the overlap 
method (see also Ref. 44). 

We used Eq. (pj to obtain a direct evaluation of the 
entropy of the M solute. The remaining difficulty, which 
is serious, is the potential for catastrophic loss of preci- 
sion in performing a subtraction to evaluate the partial 
molar excess internal energy, (d (U) /dnys) T . We 
followed- 50 * 51 ! in considering the mean binding energies 
of water molecules within a sphere surrounding the so- 
lute, and seeking saturation and plateau behavior with 
increasing radius of the sphere. We note again that the 
present calculation permitted the volume of the simu- 
lation cell to fluctuate with p (the pressure) specified, 
and we also note that the system size for the present 
calculation is only 512 molecules. We found, however, 
the mean net solvent binding energies to the solvent for 
water molecules in the second hydration shell of the so- 
lute was roughly zero to within statistical uncertainties of 
more than ten kcal/mol. In the first hydration shell that 



mean net solvent binding energy of water to the solvent 
was about one kcal/mol, positive, within large statistical 
uncertainties but decidedly less than ten kcal/mol. The 
mean unconditioned solute-solvent binding energy was 
-3.2 kcal/mol. The roughly estimated partial molar en- 
ergy is thus about -2.0 kcal/mol. This number is in good 
agreement with experiment and thus the value for the 
entropy is satisfactory as well. The first such direct en- 
tropy evaluation, for liquid waterj^ was strikingly sim- 
ple, but here the study of the temperature dependence 
of the free energy would be more efficient. An impor- 
tant physical point is that this evaluation reinforces the 
conclusion that solvent-solvent interactions in the neigh- 
borhood of a primitive hydrophobic solute can be pos- 
itive, a point emphasized recently by conclusive results 
for hard-sphere solutes.' 5 -' This point suggests limitations 
in the classic iceberg or clathrate pictures of hydropho- 
bic hydration. Of course, such values are expected to be 
strongly temperature dependent. 

FIG. H] also shows that extreme value distributions^ 
are effective in modeling the distributions -P(e), partic- 
ularly the unconditioned distribution. The high-e be- 
havior is faithfully exponential, and the Gumbel extreme 
value distribution lnP(e) = - (e+ 3.52) /0.59 + 0.52 - 
e -0+3.52)/o.59 j g accura te. This further supports the view 
that the characteristic form of the unconditioned distri- 
bution in the high-e region is due to energetic (i.e. 'ex- 
treme') interactions with a small number of molecules. 

A gamma distribution of interaction energies has been 
proposed in a different statistical thermodynamic setting 
for closing and integratin g a di fferential equation for an 
entropic equation of statePSSH A gamma distribution is 
qualitatively satisfactory here but not as accurate as the 
Gumbel distribution. As does the present work, those 
previous studies emphasized realizability in free energy 
models: free energy evaluations can benefit from restric- 
tion to legitimate probability distributions. 

A distinction from that preceding workP^I] is that we 
study here the distribution of an intensive characteris- 
tic, a possibility noted before.^ Then the central limit 
theorem is not available to force the unconditioned distri- 
bution toward a normal (gaussian) form. Alternatively, 
the present theory naturally presents a physical control 
feature, the volume of the defined inner shell, to drive 
the observed conditioned distribution toward a normal 
(gaussian) form. With minimal conditioning as here, the 
conditioned distribution can be slightly swper-gaussian. 
With aggressive conditioning, these binding energy dis- 
tributions can become SM&-gaussianPS Then an alterna- 
tive model, perhaps a beta distribution, might be more 
effective. 

Returning to the results, an important point here is 
that the chemical and fluctuation contributions depend 
weakly on r, as FIG. [3] shows. 

FIG.gshows AwmmM + ^m ( see E q s - ^ and ©)• 
The agreement between the gaussian quasi-chemical the- 
ory and the numerically exact overlap result is close ex- 
cept in the region of the dehydration barrier. We ascribe 
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FIG. 3: Notice the slight variation with radius r of the chemi- 
cal contribution kTlnpuinx =0|r) and the fluctuation contri- 
bution (5e 2 \r, n\ =Q)/2kT. The probability that one methane 
molecule of a pair has an empty inner-shell (pu(nx = 0\r)) is 
greatest just outside MM contact (see FIG.[4|, and least just 
outside the predicted dehydration barrier where intervening 
water molecules can take effective advantage of attractive in- 
teractions with each methane molecule of the pair of methane 
molecules. 



the discrepancy of the theory-overlap data comparison in 
that region to the difficulties of the present direct esti- 
mate of p(n\ = 0\r) (see Eq. ( |16[ )), particularly in view of 
the evident noise in the theoretical result. Together with 
the results of FIGS. [2] and [3j this is strong confirma- 
tion of the physical correctness of the mean-field descrip- 
tion of the effects of attractive interactions in these pm/s. 
FIGS. [3] and [4] together show that for molecular size so- 
lutes, and where the strength of attractive interactions 
with water are similar to those expressed by aliphatic 
groups, the packing contributions are the biggest single 
contribution though mean-field contribution from attrac- 
tive interactions are not negligible. The net effect here 
of attractive solute-water interactions is repulsive, as was 
guessed beforehand we expect that to be the usual case. 

Discussion of these point s, sh arpens our understanding 
of some earlier calculations 8 55 that utilized solely pack- 
ing contributions for such pair potentials of mean force. 
In a precisely denned setting, such models should be aug- 
mented to include the effects of attractive forces. Just 
as important, information theory models for co mposi te 
species that can dissociate have some peculiarities.!^ 4 ^ It 
is likely that those peculiarities combined with secondary 



but non-negligible affects of solute- water attractive inter- 
actions are part of the understanding of the striking ini- 
tial observations of agreement of the information theory 
models with simulation results. 8 




FIG. 4: The absolute solvent contribution to the potential 



of mean force between two methane molecules, see Eq. ( 14 1 



The bold line i s th e reference value obtained from the overlap 
method of Eq. ( 17 1 . The triangles are the net profile obtained 



using the gaussian-model, Eq. ( 15 1 



VI. CONCLUSIONS 

The gaussian quasi-chemical theory is physically valid 
for these pair potentials of mean force. The theory here 
is sufficiently accurate that the packing contribution may 
be inferred from simulation data on the system with 
physical interations. In the present approach, the pack- 
ing contribution does not have the logical status of a 
reference system contribution in a perturbative formula- 
tion. In particular, the chemical contribution thereby 
defined is non-perturbative. Nevertheless, packing ef- 
fects make the largest contribution, and both chemical 
and fluctuation contributions are negligible by compar- 
ison. For the methane-methane pair potential of mean 
force in water at customary conditions, mean-field ef- 
fects associated with solute-water attractive interactions 
are secondary but not negligible. Attractive solute-water 
interactions make a net repulsive contribution to these 
pair potentials of mean force. This substantiates a pre- 
vious physical guess and is consistent with intervening 
observations.^ Extreme value distributions accurately 
model the distributions observed here, which are super- 
gaussian. 
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